% This cript computes tidal power spectra for all the geomagnetic
% for the station Bangui. The idea is to study the intersting M2 power wer
% observed in this station

% Get a list of stations

S = dir('/data/backup/mnair/obs_mag_data/update_2011/secular_removed/*.nc');
a = (datenum(2011,1,1)- datenum(1995,1,1))*24*60;
t = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2010,12,31,23,59,30));
%periods = [4 4.8 6 8 11.967236 12  12.421 12.6583 23.934472 24 25.891];

periods = 12.4206012;

periods = 12.421;

x = [cos(2*pi*t/( periods(1)/24))'       sin(2*pi*t/( periods(1)/24))'];

for j = 2:length(periods),
    
    x = [x cos(2*pi*t/( periods(j)/24))'       sin(2*pi*t/( periods(j)/24))'];
    
end;


nstn = 21; % BNG Bangui station
    
ncid = netcdf.open(['/data/backup/mnair/obs_mag_data/update_2011/secular_removed/' ... 
    S(nstn).name],'NOWRITE');


X_ID = netcdf.inqVarID(ncid,'Magnetic_Field_X');
Y_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Y');
Z_ID = netcdf.inqVarID(ncid,'Magnetic_Field_Z');
z_data = netcdf.getVar(ncid,Z_ID,0,a);

latitude  = double(netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'geospatial_lat'));
longitude = double(netcdf.getAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'geospatial_lon'));

z_data = double(z_data)/10;
z_data(z_data == 99999.9 ) = NaN;
L = isnan(z_data);
k = find(L==0);
z_data(L) = 0;

netcdf.close(ncid);


%Found the the data quality for Bangui station for early years are pretty
%bad. Considering only the post 2005 data.

L = L | t' < datenum(2005,01,01);



%robust fit
%all data

spectra_all = robustfit(x(~L,:), z_data(~L));
% only nightime data
%Get the LT
t = t + (longitude/15)/24;

 L = L | (t-floor(t)  < 18/24 & t-floor(t)  > 6/24)';
 spectra_night = robustfit(x(~L,:), z_data(~L));
 

%trying out the difference between 1 component inversion and
%multip-component inversion
%Also found that robust fit is not yeilding great differences from simple
%LS fitting

% prd = 12.415;
% clear spectra_data_ls_m2_n
% for i = 1:100,
%     
% x_m2 = [ones([1,length(t)])' cos(2*pi*t/( prd/24))'  sin(2*pi*t/(prd/24))'    ... %M2
%     ];  
% spectra_night_ls_m2 = x_m2(~L,:)\z_data(~L);
% spectra_data_ls_m2_n(i) = abs(complex(spectra_night_ls_m2(2),spectra_night_ls_m2(3)));
% prd = prd + 0.0001;
% end;
% 
% prd = 12.42;
% clear spectra_data_ls_m2_n
% for i = 1:100,
%     
% x_m2 = [ones([1,length(t)])' cos(2*pi*t/( prd/24))'  sin(2*pi*t/(prd/24))'    ... %M2
%     ];  
% spectra_night_ls_m2 = robustfit(x_m2(~L,:),z_data(~L));
% spectra_data_ls_m2_n(i) = abs(complex(spectra_night_ls_m2(2),spectra_night_ls_m2(3)));
% prd = prd + 0.0001;
% end;
% 
